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THE WAVE-FRONT TRACKING ALGORITHM FOR HUGHES’ 
MODEL OF PEDESTRIAN MOTION* 

PAOLA GOATIN AND MATTHIAS MIMAULTt 

Abstract. We code the wave-front tracking algorithm for Hughes’ model of pedestrian motion 
with generalized running cost. This model displays a non-classical dynamic at the splitting point 
between the two directions of motion. The wave-front tracking scheme provides us with reference 
solutions to test numerically the convergence of classical finite volume schemes, which do not treat 
explicitly the dynamics at the turning point. 

Key words. Macroscopic pedestrian fiow models, conservation laws, eikonal equation, wave- 
front tracking, finite volume schemes. 

AMS subject classifications. 35L65 


1. Introduction. We consider a generalization of the macroscopic pedestrian 
flow model introduced by Hughes [9]. The pedestrian density p = p{t,x) is assumed 
to evolve according to the following system consisting in a scalar conservation law 
expressing mass conservation coupled with an eikonal equation giving the preferred 
direction depending on the density profile, which in the one dimensional case read: 

Pt - =0, (1.1a) 

\(l)x\=c{p), (1.1b) 

in the spatial domain O = ] —1,1[. The above equations must be completed with 
homogeneous (weak) Dirichlet boundary conditions at x = ±1, which represent the 
exits locations 


p(t,-l)=p(t,l) = 0, (1.2a) 

0(t,-l) = 0(t,l) = O, (1.2b) 

and an initial density profile p(0, •) = po ^ BV(M). 

Above, X G O is the space variable, t > 0 is the time, p = p(t, x) G [0,1] is the 
(normalized) crowd density, c = c(p) is the cost function, v{p) = {1 — p) is the mean 
velocity and we set /(p) = pv{p) = p(l — p). We assume that c : [0,1[ ^ [l,+oo[ 
is a smooth function such that c(0) = 1 and c'(p) > 0 for p G [0,1[. Recall that the 
literature usually proposes the choice 


c{p) = 


1 

v{p) 


(1.3) 


n min Ho]. 

We observe that (1.1a) can be rewritten as 


Pt - {f{p) = 0, 


(1.4) 
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and the unique viscosity solution to the Dirichlet problem (1.1b), (1.2b) is given by 
the value function of the corresponding control problem with discontinuous coefficient 
c(p), i.e. 


<p{t,x) = < 


J {pit, y)) 

[ c{p{t,y)) 

^ J X 


dy if a; < ^(t), 


(1.5) 


dy if a; > ^{t), 


where ^{t) is implicitly defined by the identity 
riit) 

c{p{t,y)) dy = 

im 


riit) cl 

/ c{p{t,y)) dy= c{p{t,y)) dy. 

J-i Jf(t) 


( 1 . 6 ) 


Therefore, equation (1.1a) can be written as a scalar conservation law with discontin¬ 
uous space-time dependent flux: 


Pt + = 0, (1.7) 

where F(t,x,p) = sgn (x — ^(t)) /(p). Observe that the position of the discontinuity 
is not a priori fixed, as assumed in previous results on conservation laws with discon¬ 
tinuous fluxes (see for example [12] and references therein), but depends non-locally 
on p itself. Therefore, known convergence results for finite volume schemes do not 
apply in this setting. 

From the mathematical point of view, problem <0) presents a non trivial cou¬ 
pling between a scalar conservation law and an eikonal equation, that poses several 
challenging questions concerning existence, uniqueness, numerical approximation and 
construction of the solutions. Some preliminary results concerning properties of so¬ 
lutions have been obtained recently in BE]. In particular, they provide a notion of 
entropy weak solutions, a discussion on their qualitative behavior, and an algorithm 
to construct locally in time a self-similar solution in a neighborhood of the turning 
point X = ^{t) (see [U Theorem 1]). 

Here we are interested in computing numerically the weak solutions of the (non- 
regularized) hyperbolic problem As remarked above, classical finite volume 

techniques for discontinuous fluxes do not apply here due to the fact that the discon¬ 
tinuity location is not a-priori fixed, but depends non-locally on the solution itself. 
Nevertheless, it appears that a direct application of classical shemes allows to recover 
the exact solutions. Since these are not known explicitly, the convergence of the finite 
volume schemes is numerically tested on accurate approximations of entropy weak so¬ 
lutions computed using the wave-front algorithm (for a detailed review of the method 
we refer to [7]). 

2. Basic properties of solutions. For sake of simplicity, we refer to the initial¬ 
boundary value problem 


pt + F{t,x,p)x = 0, 

\<i>x\ = c{p), 
p{t),x) = Po{x), 

P(i,-1) = p{t,l) =0, 

= 0, 


X €0, ,t>0, 
X €0. 


( 2 . 1 ) 


t > 0 
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Following [5], we look for the following class of solutions. 

Definition 2.1. (Entropy Weak Solutions). 

A function p G L^(0)) flBV (M+ x [0,1]) is an entropy weak solution of the 

initial-houndary value problem ( |2.1[ ) if the following Kruzkov-type entropy inequality 
holds for all k G [0,1] and all test functions G C^(M x fl;]R+); 


^-|-00 nl nl 

0 </ / {\p - k\tpt + p,k)tpx) dx dt+ / \po{x) - k\tp{0,x) dx 

Jo J-i J-i 

r-\-00 

+ sgn{k) / (/ {p{t, 1 -)) - f{k)) 1 ) dt 

Jo 

n-\-oo 

+ sgn(fc) / (/ {p{t, -1+)) - f{k)) ip{t, -1) dt 

Jo 

/• + 00 

+ 2 / f{k)f;(t,m) dt. ( 2 . 2 ) 

Jo 


Above, ^(t, X, p, k) = sgn(p — k) {F{t, x, p) — F{t, x, k)) and the boundary conditions 
are intended in weak form as in [2]. In particular, the traces of the solution at the 
boundary points must satisfy 


/ {p{t, -1+)) > f{k), for all k e [0, p{t, -1+)], 

/ {p{t^ 1-)) > /(^)> for all k e [O, p{t, 1-)]. 


This in particular implies p(t, —1+) <1/2 and p(t, 1—) < 1/2. 

Besides, by taking suitable test functions vanishing along x = ^(t), it is easy to 
recover the Rankine-Hugoniot condition 

f {p'^) + f {P~) = i{t) {P'^ - P~), (2-3) 

where we have denoted by = p^(t) = p(t,^(t)+) and p~ = p~{t) = p(t,^(t)—) 
the right and left traces of p at x = ^(t). We remark here that if p(t, •) is continuous 
at X = ^(t), then ( |2.3[ ) gives / (p (t, C(^))) = 0 and therefore p (t, f,{t)) G {0,1}. 

For further properties of entropy weak solutions, such as maximum principle, we 
refer the reader to [5]. 

3. The Riemann solver at the turning point. The argument proposed in [1] 
is based on the observation that, due to finite wave propagation speed, for every fixed 
t > 0 the solution p(t,x) is well defined away from the turning point f = ^{t) for 
t > f sufficiently small, say on [—1, ^ — J) U (^ + (5,1] for some ^ > 0. Therefore, the 
following quantity is well defined 

dt[c{p{t,y))]dy =: 

Moreover, let us assume that the the function p(f, x) is constant in a left and right 
neighbors of f, let’s say 



p{t,x) 


pL if ^ - 26 < X < ^, 
PR if ^ < X < C + 25. 


(3.1) 


Assuming T* to be constant for small times t > the solution around the turning 
point will be self-similar, behaving as a solution of a Riemann problem. This as¬ 
sumption is met when there are no interacting patterns in the whole interval [—1,1]. 
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Therefore, for small times t > t we have 

np] = ^* + \ - 1 dt[c{p{t,y))]dy. (3.2) 

[Jm Ji-s J 


The term T* depends on p(f, •) in a non-local way and may range all over M, 
independently of and pR. The solution around f will be classified according to the 
value of T*. 

Theorem 3.1. (Riemann Solver at f = C(t)). Assume that (3.1) and (3.2) 
hold, then a self-similar solution around (t, f) is defined loeally in time as follows: 

1. If pL > pR, the following eases oeeur. 

la) If 


< h^A±Ikll[cipR)+c{pL)], (3.3) 

PR - pL 

then there exists a unique intermediate value pM, with pL > Pm > 
Pr, sueh that the solution is given by the turning eurve ^ followed by a 
rarefaetion between pM ci'^d pR. 

lb) If 


fiPn) + f{PL) 


c{pr) +c(pl)] < 5'* 


pR - pL 

< -v{pl){1 + c{pl)) - v{pr){1 - c{pr)), 


(3.4) 


then a unique intermediate value pM ^ [0, Pr\ exists, sueh that the solu¬ 
tion is given by the turning eurve ^ followed by a shoek between pM cind 
pR. If the equality holds in the r.h.s. of ( |3.4| ), then pm = 0.* vaeuum 
appears between the turning point and the shoek. 

Ic) If 


-v{pl){1 + c{pl)) - v{pr){1 - c{pr)) < 

< v{pr){1 + c{pr)) + v{pl){1 - c{pl)), 

then the solution is given by a shoek of speed Pl ~ 1 followed by ^ and by 
a shoek of speed 1 — pR, the intermediate state around ^ being pM = 0. 

Id) Finally if 


T* > v{pr){1 + c{pr)) + v{pl){1 - c{pl)), 


(3.6) 


then the solution is given by a shoek between pR and pM ^ [0, Pr\ followed 
by If the equality holds in (3.6), then pM = 0, otherwise pM > 0. 

2. If pL < pRy the following eases oeeur. 

2a) If 


T* < -v{pl){1 + c(pl)) - ^(pr)(1 - c{pr)), 


(3.7) 


then a unique intermediate value pM ^ [O^Pl] exists, sueh that the solu¬ 
tion is given by f followed by a shoek between pM and pR. If the equality 
holds in (3.7), then pM = 0.* vaeuum appears between the turning point 
and the shoek. 
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2b) If 

-v{pl){1 + c{pl)) - v{pr){1 - c{pr)) < 

< v{pr){1 + c{pr)) + v{pl){1 - c{pl)), 

then the solution is given by a shoek of speed pL — ^ followed by f and 
by a shoek of speed 1 — pR, the intermediate state being pM = 0. 

2c) If 


v{pr){1 + c{pr)) + v{pl){1 - c{pl)) < 4 '* 

^ f{PR) + Upl) r . \ I / M 

< ----- Hpr) + c{Pl)], 

PR - PL 


(3.9) 


then the solution is given by a shoek between pL and pM ^ [O^Pl], fol¬ 
lowed by the turning eurve If equality holds in the l.h.s. of ( |3.9[ ); then 
Pm = 0; otherwise pM > 0. 

2d) Finally if 


pR - pL 

then there exists a unique intermediate value pM, with pR > pM > 
pL, sueh that the solution is given by a rarefaetion between pR and pM 
followed by the turning eurve 
3. If pL = pR = P*; the following eases oeeur. 

3a) If 




(3.11) 


3b) 


then a unique intermediate value pM ^ [0? P*] exists, sueh that the so¬ 
lution is given by the turning eurve f followed by a shoek between pM 
and pr. If the equality holds in ( 3.11| ); then pM = 0; vaeuum appears 
between the turning point and the shoek. 

If 


-2v{p^) < ^* < 2t;(p*), (3.12) 


then the solution is given by a shoek of speed p* — 1 followed by f and by 
a shoek of speed 1 — p*; the intermediate state aeross f being pM = 0. 
3e) Finally if 


> 2t;(p*), 


(3.13) 


then the solution is given by a shoek between pR and pM ^ [0, Pr\ followed 
by If equality holds in (3.13); then pM = 0; otherwise pM > 0. 

The proof follows closely the one of [TJ Theorem 1] and it is therefore omitted. 
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4. Numerical schemes. 

4.1. The wave-front tracking scheme. The wave-front tracking technique 
for classical scalar conservation laws consists in constructing piece-wise constant ex¬ 
act solutions of an approximated problem with piecewise constant initial datum and 
piecewise linear flux /^, coinciding with / = pv{p) on the mesh Mv = ^ 

defined by 


Af^ = (2-^Nn[0,l]) 

for z/ G N, > 0. We refer the reader to 013 for further details. We describe here 
the extension of the method to the case of the coupled problem 

Let Pq be a piecewise constant function defined by 

Po=Yl Po,j ^ 

jez 

which approximates po in the sense of the strong topology, that is 


-pollLi(n) =0’ 


and such that TV(/9o) < TV(/9o)- Above, we set xq = which is calculated from 
the cost balance equation (|1.6|) at f = 0: 


that reduces to 



c{poiy)) dy 



c{poiy)) dy, 


- Xj-i) = J2poj(^j - 

j<0 j>0 


Observe that we may have Pq q = Po i • 

For small times t > 0, a piecewise approximate solution (p^ 
structed piecing together the solutions to the Riemann problems 


to (2.1) is con- 


dtp + {sgn{x - ^o)f''{p)) =0, 
poiix < ^0, 

Pi if X > ^0, 

^{p'^ - P~) = Mp], 


p(0,x) = 


dtP + 9a, (sgn(xj 


eo)r (p)) = 0, 


p(0,x) 


pj if X < Xj, 
Pj^i if X > Xj, 


jVO. 


(4.1) 


Solving the problem by front tracking replaces the regular rarefaction fronts by 
rarefaction fans of constant values p^, j = /i,..., such that |p^ — p^-i\ =2“^, 


4.1 


separated by jump discontinuities moving with speeds Xj = 1 — pJ — Pj_i, see Fig 
Note that the solution to the Riemann problem in (E3, left, is constructed by means 


of the Riemann solver described previously in Theorem |3.1[ where the theoretic value 
of the intermediate state pM is approximated by the closest point of the mesh My. 
This introduces an error in the resolution of the linearized problem, but guarantees 
that the approximate solution keeps lying on My. The (piece-wise constant) speed 
of the approximated turning curve is then computed using the Rankine-Hugoniot 
equation (2.3). 
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Fig. 4.1. Wave-front tracking solution of with po(x) = 0 for a: G ] — 1, 0[ and po{x) = 0.9 

for X G ]0,1[; and mesh size Ap = The white curve is the trajectory x = C(0- 



Fig. 4.2. Density profile at time t = 0.8 with initial density po{x) = 0 for x G ] — 1,0[ and 
p{){x) = 0.9 for X G ]0,1[, computed using wave-front tracking scheme with mesh size Ap = 2“^®. 
The position of the turning point x = ^(t) coincides with the first jump discontinuity from the left. 
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The piecewise constant approximate solution constructed above can be pro¬ 
longed up to the fist time t > 0, where two discontinuities collide, or a discontinuity 
hits the turning curve In both cases, a new Riemann problem arises and its solu¬ 
tion, obtained in the former case with the classical Riemann solver and in the latter 
by means of Theorem |3.1[ allows to extend further in time. We point out that 

the absorbing boundary conditions (1.2a) are taken into account by simply dropping 
the waves hitting the left and right boundaries. 

This procedure will provide us with a reference solution to numerically test the 
convergence of classical finite volume schemes. Even if the theoretical convergence 
of wave-front tracking approximations to exact solutions of (13 makes the object of 
on-going research, numerical evidence of this convergence can be found in the tests 


reported in Section 5.2 


The Matlab code used for the numerical tests presented in Section can be 
downloaded at the following URL: 

http://www-sop.inria.fr/members/Paola.Goatin/wft.html, 

4.2. Finite volume schemes. In this section, we describe the algorithm used 


in [5] for numerical simulations of the model (1.1). Finite volume schemes are com¬ 
monly used to compute numerically solutions of the 2D version of the model, see for 
example um and references therein for more details and results. 

Here we aim at investigating the behavior of classical schemes near the turning 
curve, and whether they manage to capture correctly the behavior of the solution 
constructed via the wave-front tracking algorithm. 

Given an initial datum p(0,x) = po{x) and homogeneous Dirichlet boundary 
conditions (1.2), we solve ([H in an iterative manner at each time step, i.e. 

1. Given p, solve the eikonal equation ( |l.lb ) by the fast sweeping method, see 


Section 14.2.11 


2. Given 0, solve the non-linear conservation law (1.1a) using Godunov or Ru¬ 
sanov scheme. 

The domain [—1,1] is divided into N uniform cells Ij = [:rj_i /25 ^j+ 1 / 2 ] with centers 
at points Xj = jAx, with Ax = 2/N. The explicit algorithm used to generate the 
approximations p^ was introduced by Towers in m and is written in conservation 
form 


pr^=p] - 


At^ 

Ax 




j 2 


,)■ 


(4.2) 


where Af^ = is chosen to satisfy the following GLF condition 

Ax 


Ae < 0.5 


max|maxj \f'{pp\, |C"|} 


In the above formula, the coefficient 0.5 is chosen to avoid interactions of ^ with 
the cell boundaries, and |^| is estimated at each time step by deriving the implicit 
expression (1.6), which gives 


m 


r^it) p1 

{c{p~) + c{p+)) =- c{p{t,y))t dy+ c{p{t,y))t dy. 

J-i Jm 


From the above identity we recover the (sharp) upper bound 


in < 2 
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In (4.2), we set ^^± 1/2 


sgn { 4 >l{xj±i/ 2 ))-, where 


<t>x{Xj+l/2) 


Ax 


The numerical flux is chosen to be monotone and consistent, 

i.e. h{p, p) = f{p) = pv{p). In order to maintain the monotonicity of the scheme, we 
transpose the arguments when %+i /2 changes sign, i.e. 


I h{pj,pj+i) iffc^.+ i>0, 

hi+l = { (4-3) 

[ h{pj+i,pj) if fcj+i < 0. 

4.2.1. The Fast Sweeping algorithm jl5j . 

Discretization. The method uses a Godunov upwind difference scheme to discretize 
the partial differential equation at interior points of the domain: 

= c{pi)‘^h^, j = 2 ,N -1, (4.4) 

where (j)'} ~ and and 


(:r)+ = 


X, X > 0, 

0, X < 0. 


Initialization. We assign exact values at boundary grid points x = ±1 to enforce 
boundary condition 0(±1) = 0. We assign sufficiently large positive values at all other 
grid points. These points will be updated later. 


Iterations. At each grid point Xj whose value is not fixed during the initialization. 


compute the solution of (4.4), denoted by 0, from the current values of its neighbors 
and then update (pjto be the smaller one between 0 and 0^, i.e. = 

min((/)9^‘^, 0). We sweep the whole domain with two alternating orderings repeatedly: 


(l)j = l:7V, (2)j = 7V:l, 


and we stop when the fixed error threshold is met. 

4.2.2. Godunov scheme. The Godunov scheme [6] for a conservation law of 
type Ut + F{u)x = 0 is classically obtained by using the exact solution with piecewise 
constant initial data. The numerical flux is h{u,v) = F where 7l{0;u,v) 

is the solution of the Riemann problem with left and right states u and v evaluated 
at X = 0, where a jump on the initial data occurs. 

The Godunov flux is given by 


(4.5) 


We remark that in the case of problem an explicit Riemann solver at the sign 


change interface would require to code all the possible cases detailed in Theorem 3.1 


and would be very expensive. Aiming to have just a qualitative representation of the 
solution, in our computations we force to be at the middle of the cell it belongs to 
at each time step, i.e. = Xj, for j such that sgn{(l)x{xj- 1 / 2 )) > «5^^(0cc(^j+i/2))- 
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4.2.3. Rusanov scheme. In [14], the author proves convergence of Engquist- 
Oscher scheme for conservation laws with discontinuous flux. These results can be 
used in the second step to solve the conservation law. For our simulations we have 
used a variant introduced by Rusanov whose flux is given by 

hiu,v) = ^{f{u)+f{v)) + ^max{\f{u)\,\f{v)\}{u-v), (4.6) 

and is known to be robust. 


5. Convergence results. In this section we provide some numerical tests show¬ 
ing the convergence of the wave-front tracking and finite volumes schemes described 
in the previous Sections |4.1[ |4.2[ To this end, we choose an initial datum po that 
gives a non-trivial behavior at the turning curve. In fact, if the solution p of (2.1) is 
continuous at x = ^(t), that is, p(t,^(t)—) = 0 = p(t,^(t)+), we deal with classical 
solutions on each side of the turning curve, and the convergence is standard. 

In the following, we consider the Riemann-type initial datum 


Po{x) 


0 if — 1 < X < 0, 
0.9 if 0 < X < 1, 


(5.1) 


and the cost function c(p) = ^-/v{p) = 1/(1 — p). We run simulations up to time 


T = 3, so that the whole mass lefts the domain Vt and p{T^x) = 0, see Fig. 4.1 


5.1. Data processing. Wave-front tracking solutions are not defined on a stan¬ 
dard cartesian grid in space and time, but rather as constant values on polygonal 


regions delimited by wave fronts (see Fig. 5.1). In order to compare the approximate 
solutions obtained via wave-front tracking at different mesh sizes, and finite volume 
approximations with the wave-front tracking reference solution, we have to convert 
these data sets to values defined on a reference cartesian grid. 

In what followsjve denote by Xwft the data obtained from the wave-front track¬ 
ing procedure, by Xwft the processed data and by Xpy those obtained from finite 
volumes schemes. Values for Pwft obtaine d by taking the corresponding values 
of Pwft ^t the cartesian grid nodes, see Fig. 5.1 Finally, we operate a L^-norm 
comparison, according to the following formula : 


T 1 N J 

Err{p,Ax)= [ [ = V V KpvVft)" “ (PFy)”| ArcAi", 

Jo J-l n=li=l 


where v refers to the wave-front tracking mesh size and Xx denotes the space-mesh 
size of the cartesian grid used for the finite volume scheme. 

We remark that the above procedure introduces an approximation error in the 
computation of the L^-error, since 



\PWFT- PFv\dxdt ^ ! L \pWFT ~ pFv\dxdt. 


5.2. Convergence of Wave Front Tracking. The theoretical proof of exis¬ 
tence and uniqueness of solutions of (2.1) by means of the wave-front tracking method 
is currently under investigation. Meanwhile, we can produce numerical evidence of the 
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Fig. 5.1. Definition of pw ft cit the point {tpv t^fv) of the referenee eartesian grid: after 
finding the last interaction time t'^ < typ in PWFT, one has to compare the positions of the 
corresponding interfaces Xi{tFv) = -\-Si{tFV —t^) with respect to xpv, where Si = 1 — pi — piyi 
are the speeds of the wave-fronts. Here pwFT = pff = • 


convergence of WFT scheme by showing that the -error between two subsequent 
discretization meshes 2“^ and 


— IIpivft ~ Pwft\ 


Li 


decreases to zero as u increases. Table ] 5.1 [ and Fig. 
is a Cauchy sequence in L^. 


5.2 


show that the sequence {Pwft} 


u 

Ap 


5 

2-^ 

4.280e - 2 

6 

2-^ 

2.164e-2 

7 

2-^ 

6.141e-3 

8 

2-^ 

5.048e - 3 

9 

2-9 

1.755e-3 

10 

2-10 

2.091e-3 

11 

2-11 

4.305e - 4 

12 

2-12 

4.347e - 4 


Table 5.1 

L1-error Cjy for wave-front tracking method between two subsequent discretization meshes 2~^ 
and The comparison is done on a cartesian grid with Ax = 10“^ and At = 0.5Aa:. 


5.3. Convergence of finite volume schemes. Once assessed the convergence 
of the wave-front tracking scheme, we take as reference “exact” solution the one 
obtained with density mesh size 5 = 2“^ = , and we compare it to the approxi¬ 

mations computed by Godunov and Rusanov schemes for diff eren t space mesh sizes, 
computing the L^-norm of the difference. The following Table 5.1 and Fig. 5.3 report 
the values of the -error computed on the time interval [0, T] with T = 1.2, in order 
to focus on the non-classical behavior of the solution. Indeed, at T = 1.2 vacuum has 
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Fig. 5.2. -error between two subsequent di seret ization meshes 2 ^ and 2 ^ ^, for u = 
5,..., 12, eorresponding to the initial datum given by 


Ax 

Errc 

In(^rrG')/ In(Ax) 

Erru 

ln{ErrR)/ ln{Ax) 

1/50 

7.24e - 2 

-0.66 

7.44e - 2 

-0.67 

1/100 

4.56e - 2 

-0.66 

4.68e - 2 

-0.67 

1/250 

2.49e - 2 

-0.66 

2.55e - 2 

-0.67 

1/500 

1.52e-2 

-0.67 

1.55e-2 

-0.67 

1/1000 

9.03e - 3 

-0.68 

9.12e-2 

-0.68 

1/1500 

6.66e — 3 

-0.69 

6.62e - 3 

-0.68 


Table 5.2 


L^-norm of the error for Godunov and Rusanov sehemes for depending on the spaee step 

Ax, and eorresponding eonvergenee order. Above, Ewq stands for the error in Godunov seheme, 
ErvR indieates the error for Rusanov flux. 


appeared around the turning point and the solution has become classical. For a 
set of space meshes ranging from Ax = 1/50 to Ax = 1/1500, we see a clear linear 
decreasing in logarithmic scale, showing numerically a convergence order of about 
0.67 for both Godunov and Rusanov schemes, without an explicit treatment of the 
solution’s behavior at x = ^(t). Fig. 5.4 shows the profiles of the three numerical 


approximations in a neighborhood of the turning point. We observe that Godunov 
scheme displays small oscillations close to x = ^{t). 


5.4. More general initial data. We provide here the results of simulations 
performed considering a more general initial datum. We have taken 


^ 0.8 if -0.8 < T < -0.5, 
0.6 if — 0.3 < X < 0.3, 
0.9 if 0.4 < X < 0.75, 

0 elsewhere. 


Po = < 


(5.2) 
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Fig. 5.3. 'L^-norm of the error in logarithmie seale for mesh sizes Ax = 1/ 50,..., 1/1500; for 
Godunov and Rusanov schemes corresponding to the initial datum given by (|5.1|>. 



Fig. 5.4. Zoom o f the numerical approximations in a neighborhood of the turning point x = ^(t) 
for the initial datum at time t = 0.8. The wave-front tracking profile is given by the mesh 

Ap = 2“^®; and the finite volume space step is Ax = 1/1500. 


as in [H Figure 6] in order to simulate the behavior of three mixed-density groups. 
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Goatin and Mimault 





Fig. 5.5. 
po given by 
Aa; = 10-3. 


Wave-front traeking (top), Godunov (middle) and Rusanov (bottom) sehemes for 
I The wave-front traeking mesh is Ap = and the finite volume spaee step is 


We can see that the two finite volume simulations are evolving like the wave- 
front tracking one, modeling two turning phenomena. We can also observe that the 
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diffusion of Rusanov counterbalances severely Godunov’s oscillations. 

6. Conclusions. A crowd evacuating a corridor can be modeled by a system of 
partial differential equations in one space-dimension, coupling a scalar conservation 
law describing conservation of mass and an eikonal equation assigning the direction 
of motion depending on the density profile. The system can be rewritten as a scalar 
conservation law with space-discontinuous flux function, for which the discontinuity 
location depends non-locally on the density profile. 

Since exact entropy weak solutions cannot be computed theoretically for general 
initial data, we have coded an adapted wave-front tracking scheme, which is reason¬ 
ably assumed to give a good approximation of problem’s solutions. These wave-front 
tracking approximations are used as reference solutions to show numerically the con¬ 
vergence of classical finite volume schemes, which do not treat explicitly the dynamic 
at the turning point. 

These results offer a sound basis to the use of classical finite volumes schemes 
for the computation of solutions of Hughes’ model in two space dimensions, see for 

example mm- 
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